##Copy the part starting with source() into the console and run it there - this will create a log file that includes the output

#source("~/work/November_2024/November_2024_PanelMatchScript_Robustness_CoreNA.R", echo = T, max.deparse.length=10000)
rm(list=ls(all=TRUE))
sink("~/work/November_2024/R_out/log_files/full_log_CoreNA_ROBUST.log", split = T)


# ======================================================================
# Project:    Closed borders, closed minds? COVID-related border closures,
#             EU support and hostility towards immigrants
#
# Script:     Appendic Figure C1 + Figure C2
#
# Authors:    Lisa Herbig (l.j.herbig@uva.nl)
#             Asli Unan (a.unan@uva.nl)
#
# Date:       2025-03-24
# ======================================================================

# Description:
# This script generates Appendix Figure C1 and Figure C2, illustrating the average treatment effect on attitudes towards refugees and attachment with Europe with core regions being NA. 

# ----------------------------------------------------------------------
# Load libraries
# ----------------------------------------------------------------------
library("did2s")
library("dplyr")
library("foreign")
library("ggplot2")
library("haven")
library("lmtest")
library("lubridate")
library("modelsummary")
library("PanelMatch")
library("panelView")
library("readr")    
library("tidyr")
library("tidyverse")


## Read the files in
df_iso_work_dd = readRDS("~/work/November_2024/R_out/df_iso_work_dd_25.rds")

# Appendix Figure C1 -----------------------------------------------------------

###CoreNA - Refugee
output_folder = "~/work/November_2024/R_out/robustness_plots/"

PM.results <- PanelMatch(lag = 2, time.id = "weekyear_i",
                         unit.id = "iso2.y_i", 
                         treatment = "any_closure_binary", 
                         refinement.method = "mahalanobis", 
                         data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                         outcome.var = "immigration_pca_n",
                         match.missing = TRUE,
                         matching = TRUE,
                         covs.formula = ~ East + age + HOMEOFFICE + female + unemployed + working_abroad + university + EU_residence  + 
                           kr_foreigner + CONTACTABROAD + Voting_2018_LeftRight+ mode_phone + iso2_pop+  I(lag(kr_inz_rate, 1:4)),
                         lead = 1:4, forbid.treatment.reversal = FALSE)

PE.results1 <- PanelEstimate(sets = PM.results, data = as.data.frame(df_iso_work_dd), se.method = "conditional", pooled = FALSE)
pdf(paste0(output_folder, "SATDEM_reverse_any_closure_binary_Placebo_ROBUST_25", ".pdf"))

plot(PE.results1)

summary(PE.results1)

# Appendix Figure C1 -----------------------------------------------------------

###CoreNA -EU
output_folder = "~/work/November_2024/R_out/robustness_plots/"

PM.results <- PanelMatch(lag = 2, time.id = "weekyear_i",
                         unit.id = "iso2.y_i", 
                         treatment = "any_closure_binary", 
                         refinement.method = "mahalanobis", 
                         data = as.data.frame(df_iso_work_dd), qoi = "att" ,
                         outcome.var = "con_europe_n",
                         match.missing = TRUE,
                         matching = TRUE,
                         covs.formula = ~ East + age + HOMEOFFICE + female + unemployed + working_abroad + university + EU_residence  + 
                           kr_foreigner + CONTACTABROAD + Voting_2018_LeftRight+ mode_phone + iso2_pop+  I(lag(kr_inz_rate, 1:4)),
                         lead = 1:4, forbid.treatment.reversal = FALSE)

PE.results1 <- PanelEstimate(sets = PM.results, data = as.data.frame(df_iso_work_dd), se.method = "conditional", pooled = FALSE)
pdf(paste0(output_folder, "Con_EU_reverse_any_closure_binary_CoreNA_ROBUST_25", ".pdf"))

plot(PE.results1)

summary(PE.results1)


sink()